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Abstract 

In  this  paper,  we  present  an  extension  of  goal-oriented  error  estimation  and  adap¬ 
tation  to  the  simulation  of  multi-scale  problems  of  molecular  statics.  Computable 
error  estimates  for  the  quasicontinuum  method  are  developed  with  respect  to  spe¬ 
cific  quantities  of  interest  and  an  adaptive  strategy  based  upon  these  estimates  is 
proposed  for  error  control.  The  theoretical  results  are  illustrated  on  a  nanoindenta¬ 
tion  problem  in  which  the  quantity  of  interest  is  the  force  acting  on  the  indenter. 
The  promising  capability  of  such  error  estimates  and  adaptive  procedure  for  the 
solution  of  multi-scale  problems  is  demonstrated  on  numerical  examples. 

Key  words:  Molecular  statics,  goal-oriented  adaptive  modeling,  error  estimation, 
multi-scale  problems,  quasicontinuum  method. 


1  Introduction 

Computational  methods  for  the  study  of  multi-scale  phenomena  have  become 
a  prominent  area  of  research  in  computational  science.  Indeed,  computing  ca¬ 
pabilities  have  reached  a  point  where  atomistic  simulations  using  quantum 
mechanical,  atomistic  potential,  and  mesoscopic  and  continuum  models  can 
be  coupled  concurrently  to  study  physical  problems  of  an  inherent  multi-scale 
nature  [2,8].  Development  of  such  methods  is  of  particular  interest  for  the 
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study  of  the  mechanics  of  materials  including  fracture  phenomena,  nanoin¬ 
dentation,  atomic  friction,  etc.,  to  name  just  a  few  [19,1].  However,  many  of 
the  existing  methods  to  date,  to  the  best  knowledge  of  the  authors,  lack  some 
analysis  of  the  error  incurred  by  coupling  such  models.  Furthermore,  conver¬ 
gence  analysis  and  comparison  studies  of  these  methods  seem  to  be  scarce  and 
not  fully  addressed  in  the  literature.  In  this  paper,  we  present  an  application 
of  goal-oriented  error  estimation  and  adaptive  modeling  to  a  model  nanoin¬ 
dentation  problem  to  partially  address  some  of  these  issues.  These  ideas  draw 
upon  work  in  [11]  where  error  estimates  for  quantities  of  interest  are  derived. 
Ideas  of  goal-oriented  adaptive  modeling  come  from  [13]  and  references  therein. 
This  goal-oriented  modeling  methodology  has  been  successfully  applied  to  the 
study  of  heterogeneous  elastostatics  and  elastodynamics,  random  heteroge¬ 
neous  materials,  as  well  as  the  study  of  linear  lattice  models  [14,16,12],  See 
also  [13]. 

In  this  study,  an  atomistic  model  based  upon  potentials  of  the  embedded- 
atom  method  (EAM)  is  used  as  a  base  model  to  simulate  the  nanoindentation 
of  a  thin  aluminum  him  [4,5].  The  target  problem  was  also  studied  in  [19]. 
Surrogate  models  are  generated  using  the  quasicontinuum  method  (QCM) 
[20,18].  Error  estimates  in  a  quantity  of  interest  are  derived  and  an  adaptive 
modeling  scheme  is  implemented  in  the  freely  available  QCM  code  [10].  A  brief 
summary  of  some  of  our  results  given  in  this  paper  were  reported  in  the  survey 
article  [13].  Here  we  give  full  details  of  an  analysis  of  multi-scale  modeling  in 
which  the  coarse-scale  modeling  is  implemented  using  the  QCM. 

The  paper  is  organized  as  follows:  following  the  introduction,  we  present  in 
Section  2  the  base  model  for  molecular  statics  problems,  derive  a  surrogate 
model  based  on  the  quasicontinuum  method,  and  describe  a  practical  example 
that  deals  with  the  nanoindentation  of  a  thin  him  aluminum  crystal.  Section  3 
is  devoted  to  the  derivation  of  error  estimates  with  respect  to  quantities  of  in¬ 
terest.  These  estimates  approximate  the  modeling  error  between  solutions  of 
the  base  and  surrogate  models.  In  Section  4,  an  adaptive  strategy  is  proposed 
for  the  control  of  the  modeling  error  by  subsequent  enrichment  of  the  surrogate 
model.  Performance  of  the  error  estimator  and  adaptive  strategy  is  demon¬ 
strated  on  the  nanoindentation  problem  described  in  Section  2.  We  finally  give 
some  concluding  remarks  in  Section  5. 

2  Molecular  statics  model 

In  this  section,  we  consider  the  problem  of  determining  static  equilibrium 
configurations  of  a  regular  lattice  of  N  atoms.  The  base  problem  is  obtained 
by  minimizing  the  potential  energy  of  the  system  consisting  of  all  atoms  in 
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the  lattice.  In  many  applications,  N  can  be  very  large  and  the  base  problem 
is  often  intractable.  In  order  to  reduce  its  complexity,  we  consider  here  the 
use  of  an  approximation  method  such  as  the  quasicontinuum  method  (QCM) 
[18-20].  In  recent  years,  QCM  has  become  a  popular  approach  for  constructing 
surrogate  problems  that  retain  only  a  small  number  of  active  atoms  during 
the  simulations.  In  that  sense,  QCM  can  be  viewed  as  a  model  reduction 
procedure. 

2.1  The  base  problem.  Let  C  be  a  regular  lattice  of  N  atoms  in  d  —  2 
or  3.  The  positions  of  the  atoms  are  given  in  the  reference  configuration  by  the 
vectors  xt  G  Rd,  i  —  1, . . . ,  N.  When  the  lattice  is  subjected  to  a  deformation 
<p  :M.d  ^  M.d,  the  atoms  move  to  the  new  positions 

Xi  =  I p(xi )  =  Xi  +  Ui,  i  =  1, . . . ,  N  (1) 

where  Ui  is  the  displacement  of  atom  i.  We  assume  that  the  lattice  in  the 
reference  configuration  covers  the  region  fl,  where  Id  is  an  open  bounded  set 
of  Wl  with  boundary  di 2.  We  also  assume  that  the  atoms  lying  on  d£2  are  all 
ascribed  essential  boundary  conditions  in  the  form 

^  =  gt,  Wxi  G  dn  (2) 

with  gL  G  Wl.  Other  boundary  conditions  will  be  considered  in  the  nanoin¬ 
dentation  application.  We  deliberately  choose  to  restrict  ourselves  to  this  case 
in  the  presentation  of  the  theoretical  results  as,  otherwise,  it  would  make  the 
exposition  rather  cumbersome  without  adding  to  the  understanding  of  the 
methodology.  Let  Na  be  the  number  of  atoms  inside  the  domain  and  Nk  the 
number  of  atoms  on  dfl  such  that  N  =  Na  +  Nf}.  Henceforth,  we  shall  use  the 
convention  that  the  interior  atoms  be  numbered  from  1  to  Na  and  the  bound¬ 
ary  atoms  from  Na  +  1  to  N.  We  will  consider  the  finite-dimensional  vector 
spaces  V  =  (SLd)N  and  Vq  =  ( Wl)Na .  In  what  follows,  we  will  conveniently 
use  the  notation  u  =  (iti,  u2 , . . . ,  %),  u  G  V,  to  refer  to  the  displacements 
of  the  collection  of  N  atoms.  Similarly,  u  G  Vo  is  the  set  of  displacements 
u=  {u1,u2,..HuN  J. 

Let  a  state  of  the  system  of  N  atoms  be  described  by  the  displacements  u  G  V. 
The  total  potential  energy  of  the  system  is  assumed  to  take  the  form 

N  N 

e(«)  =  -E/rt  +  EW  (3) 

i=  1  k=l 

where  ft  is  the  external  load  applied  to  an  atom  %  and  Ek{u)  is  the  energy 
of  atom  k  determined  from  inter-atomic  potentials.  Explicit  description  of  Ek 
will  be  given  below. 
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The  goal  of  molecular  statics  is  to  find  the  equilibrium  state  u  E  V  that 
minimizes  the  total  potential  energy  of  the  system,  i.e. 

E{u)  =  inf  E(v)  (4) 

v£V 

Vi  =  gi  on  dQ 


The  constrained  minimization  problem  is  then  equivalent  to  finding  u  E  V 
such  that: 


N 


E 


fi , 


Ui  =  gu 


i  =  l,...,Na 
i  =  Na  +  1, . . . ,  N 


(5) 

(6) 


where  d/dui  is  the  gradient  vector  with  respect  to  each  component  ui ^  l  = 
1, . . . ,  d  of  the  displacement  vector  tq,  i.e.  d/dui  =  ( d/du\ ^ . . . ,  d/dudy). 


A  variational  formulation  of  the  above  problem  is  obtained  by  multiplying  the 
Na  equations  in  (5)  by  arbitrary  vectors  v  e  Vq  so  that  the  problem  reads 


Find  u  E  V  such  that 


B(u ;  v)  =  F(v), 

Ui  =  gi, 


Vv  E  V0 

i  =  Na  +  1, . . . ,  N 


(7) 


where  the  semilinear  form  and  linear  form  F(-)  are  defined  for  any 

u  E  V  and  v  E  V0  as 


Na 


B(u;v )  =  ^ 


N 

E 


d  Ek 


du 

=i  U=i  ULLi 


u 


V,- 


Na 


n»)  =  EA- 


Vi 


i=  1 


(8) 


Note  that  Problem  (7)  is  nonlinear  in  u  and  linear  in  v.  We  assume  that  there 
exist  solutions  and  that  it  can  be  solved  by  a  quasi-Newton  method  (see  [17] 
for  details). 


2.2  The  surrogate  problem  by  the  quasicontinuum  method.  We  describe  in 
this  section  the  main  features  of  QCM.  The  reader  is  referred  to  [18,17,9] 
for  a  detailed  exposition.  The  objectives  of  the  method  can  be  summarized  as 
follows:  (i)  to  dramatically  reduce  the  number  of  degrees  of  freedom  from  Nxd, 
and  (ii)  to  substantially  reduce  the  cost  in  the  calculation  of  the  potential 
energy  by  computing  energies  only  at  selected  sites.  In  addition,  the  use  of 


4 


adaptive  approaches  for  automatic  selection  of  the  degrees  of  freedom  can 
allow  QCM  to  capture  the  critical  deformations  of  the  lattice  in  an  efficient 
manner. 

The  initial  step  of  the  method  consists  in  choosing  a  set  of  R  <C  N  represen¬ 
tative  atoms,  the  so-called  “repatoms”,  and  in  approximating  u  G  V  by  the 
reduced  vector  uq  e  W  =  (Rd)R.  The  displacements  uq  represent  the  active 
degrees  of  freedom  of  the  system  and  the  repatoms  are  conveniently  identified 
with  the  nodes  of  a  finite  element  triangulation  Vh  of  Q.  The  displacements  of 
the  (N  —  R )  “slave”  atoms  are  then  interpolated  from  u0  by  piecewise  linear 
polynomials  defined  on  the  triangular  mesh.  Let  (j>r,  r  =  1, R,  denote  the 
basis  functions  (the  hat  functions)  associated  with  Vh  and  let  Uh  be  the  finite 
element  vector  function  such  that 

R 

Uh(x)  =  ^2  u0^r(j)r(x) ,  Vx  e  Q  (9) 

r=  1 

Uo  =  (u0ti,  u0,2,  ■  ■  ■ ,  uo,r)  £  W.  The  displacements  of  the  N  atoms  in  the 
lattice  can  clearly  be  evaluated  from  uq  as 

u°  =  uh(xi),  i  =  (10) 

so  that  u°  =  . . .  ,u°N)  G  V.  This  extension  operator  will  be  referred 

to  as  77  :  W  — »  V  such  that  tv u0  =  u°.  In  a  similar  manner,  dehning  Ra  and 
Rb  as  the  number  of  repatoms  lying  in  the  interior  of  the  lattice  and  on  the 
boundary  <90,  respectively,  and  letting  ILq  =  ( Rd)Rai  we  also  introduce  the 
extension  operator  7r0  :  Wo  — >  Vo-  The  reduced  vector  t vu0  could  be  used  to 
approximate  the  total  potential  energy 

N  N 

E(tvuo)  =  -  {iruo)i  +  J2  Ek(nu o)  (11) 

i= 1  k= 1 

but  such  a  calculation  would  still  be  very  prohibitive  as  all  N  atoms  need  to 
be  visited  in  order  to  sum  up  the  atomic  site  energies. 

The  second  step  of  the  QCM  is  thus  concerned  with  and  efficient  scheme  to 
approximate  the  total  energy  E(tvUq).  The  main  motivation  here  is  to  estimate 
the  potential  energy  by  summing  only  over  the  repatoms  such  that 

R  R 

E(7VUo )  «  E0(u0 )  =  -  nr  f  0,r  '  u0,r  +  J2  nrEr(u0 )  (12) 

r= 1  r=  1 

where  nr  is  an  appropriate  weight  function  associated  with  repatom  r  so  as  to 
account  for  all  atoms  in  the  lattice,  i.e.  J2r  nr  =  N,  and  fQr  is  the  averaged 
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external  force  acting  on  repatom  r.  In  the  QCM,  the  calculation  of  the  energies 
nrEr(u0 )  is  done  in  one  of  two  ways,  depending  upon  whether  a  repatom  is 
considered  either  “local”  or  “nonlocal” .  The  attribute  “local”  refers  here  to  the 
fact  that  the  energy  at  a  point  in  the  continuum  depends  on  the  deformation 
at  that  point  only  and  not  on  its  surroundings.  Let  Ric  denote  the  number 
of  local  repatoms  and  Rni  the  number  of  nonlocal  repatoms,  R  =  Ric  +  Rni 
The  atomistic  energies  are  now  separated  into  local  and  nonlocal  contributions 
such  as: 

R  Ric  Rnl 

E  nrEr(u0)  =  E  nrElr°c(u0)  +  £  nsE?(u0)  (13) 

r=  1  r= 1  s=l 

Note  that  if  Rni  =  0,  the  method  is  called  the  local  QCM ,  and  if  Ric  = 

0,  the  nonlocal  QCM.  Otherwise,  the  method  is  referred  to  as  the  coupled 
local/nonlocal  QCM.  We  shall  only  consider  the  latter  in  what  follows. 

Local  formulation:  The  local  formulation  makes  use  of  the  Cauchy-Born  rule  [7] 
to  compute  the  sites  energies.  The  Cauchy-Born  rule  postulates  that  when  a 
crystal  is  subjected  to  a  small  linear  displacement  of  its  boundary,  then  all 
interior  atoms  are  deformed  following  this  displacement.  In  particular,  this 
means  that  every  atom  in  a  region  experiencing  a  uniform  deformation  gradi¬ 
ent  has  the  same  energy.  Since  the  QCM  uses  piecewise  linear  finite  elements, 
the  deformation  gradient  is  uniform  within  each  element  and  the  energy  in  an 
element  can  be  calculated  by  computing  the  energy  of  one  atom  only  in  the 
deformed  state.  Then  the  energy  El°c  is  given  by: 

K  I< 

nrElr°c(u0 )  =  E  K  S(Fe),  nr  =  E  K  (14) 

e=l  e=l 

where  S(Fe )  is  the  energy  of  a  single  atom  under  the  deformation  gradient 
Fe.  Here  K  is  the  number  of  elements  surrounding  the  repatom  r  and  ner  is 
the  number  of  atoms  from  nr  that  actually  live  in  element  e. 

Nonlocal  formulation:  In  this  formulation,  the  energy  is  accurately  approx¬ 
imated  by  explicitly  computing  the  energy  of  the  nonlocal  repatoms,  i.e. 
E'flfaf)  =  E8{uq).  In  other  words,  if  Ric  =  0  and  in  the  limit  case  where 
every  atom  in  the  lattice  is  made  a  repatom,  that  is  Rni  =  N,  then  nr  =  1, 
r  —  1, ...  ,N,  and  the  problem  becomes  equivalent  to  the  base  problem. 

Remark  1  ( Ghost  forces)  The  coupling  of  nonlocal  and  local  representative 
atoms  leads  to  spurious  forces,  so-called  “ghost  forces” ,  near  interfaces  of  local 
and  nonlocal  repatoms.  The  issue  is  that  the  energy  calculated  at  a  nonlocal 
repatom  may  be  influenced  by  the  displacement  of  a  local  repatom  nearby,  while 
the  converse  may  not  be  true.  Therefore  the  approximation  of  the  energy  by 
the  coupled  local/nonlocal  approach  yields  non-physical  forces  at  the  interface 
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of  the  local  and  nonlocal  regions.  A  solution  to  this  issue  has  been  devised  by 
adding  corrective  forces  to  balance  ghost  forces  ( see  e.g.  [17]). 

Remark  2  (Local/nonlocal  criterion)  The  selection  of  representative  at¬ 
oms  as  local  or  nonlocal  is  based  upon  the  variation  of  the  deformation  gradient 
on  the  atomic  scale  in  the  vicinity  of  the  atoms.  A  repatom  is  made  local  if  the 
deformation  is  almost  uniform,  nonlocal  if  the  deformation  gradient  is  large. 
In  the  QCM,  the  deformation  gradients  are  compared  element  to  element  by 
computing  the  differences  between  the  eigenvalues  of  the  right  stretch  tensor 
u  =  Vftf  in  each  element,  F  —  V</>  being  the  deformation  gradient  in  the 
elements. 

Using  the  approximation  (12)  for  the  potential  energy,  the  minimization  prob¬ 
lem  for  the  QCM  consists  of  finding  u0  e  W  such  that 

E0(u0)  =  min  E0(v)  (15) 

v  G  W 

vi  —  9i  on  <9^2 


This  can  be  rewritten  in  variational  form  as 


Find  u0  G  W  such  that 

B0(u0]v)  =  F0(v), 

«0,i  =  9ii 


Vu  G  W0 

i  =  Ra  +  1,  •  •  • ,  R 


(16) 


where  the  semilinear  form  U0  ( • ;  • )  and  linear  form  F0  ( • )  are  now  given  by 


Ra 


Bq{u]v)  =  Y) 


i= 1 


d  Er 


u) 


Vi 


Ra 


Fo(v)  =  Y  ni  fo.i  ■  vi 


1=1 


(17) 


Remark  3  (Adaptivity)  In  [9],  it  is  proposed  that  an  automatic  mesh  adap¬ 
tion  technique  be  used  to  add  and  remove  representative  atoms  “on  the  fly”, 
in  order  to  capture  the  fine  features  during  the  simulation.  The  criterion  for 
adaptivity  is  based  upon  the  derivation  of  an  error  indicator  similar  to  that  of 
Zienkiewicz  and  Zhu  [21]  for  the  finite  element  method.  The  error  indicator  is 
calculated  over  each  element  h2e  as 


■'te'  =  i/rn rr  /  (f-f)-.(f-  F)dx 

V  |  “if  |  JnK 

where  Q /^  |  is  the  volume  of  element  K ,  F(uq)  the  deformation  gradient  ob¬ 
tained  from  the  QC  solution  u0  ( piecewise  constant),  mid  F  is  a  recovered 
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smooth  deformation  gradient  obtained  by  a  L2-projection  of  F(u0)  onto  the 
finite  element  space  (piecewise  linear).  In  the  adaptive  strategy,  the  elements 
that  exceed  a  prescribed  tolerance  are  marked  for  refinement.  This  adaptive 
strategy  will  be  used  as  is  to  compute  the  solution  u0  and  the  “overkill”  so¬ 
lution  of  the  problem.  Our  goal  here  is  to  propose  an  alternative  method  that 
automatically  adapts  the  solution  process  by  controlling  errors  in  quantities  of 
interest. 

2.3  An  application  example.  Numerical  simulations  to  illustrate  the  perfor¬ 
mance  of  the  error  estimator  and  adaptive  strategy  will  be  performed  on  the 
nanoindentation  problem  suggested  by  Tadmor  et  al.  [15,17,19].  This  exam¬ 
ple  is  actually  provided  as  a  model  example  accompanying  the  open  source 
software  package  [10]. 

In  this  example,  a  thin  him  of  aluminum  crystal  is  indented  by  a  rigid  rectan¬ 
gular  indenter,  infinite  in  the  out-of-plane  direction,  as  depicted  in  Fig.  1.  The 
dimensions  for  the  block  of  crystal  are  2000  x  1000  A2  (A=Angstrom)  in  the 
[111]  and  [110]  directions  of  the  crystal.  The  crystal  rests  on  a  rigid  support 
so  that  homogeneous  boundary  conditions  Ui  =  0  are  prescribed  for  those 
atoms  i  located  at  y  =  0.  The  remaining  boundary  conditions  are  enforced  as 
follows:  at  x  —  0  and  x  =  2000,  homogeneous  essential  boundary  conditions 
are  prescribed  in  the  x-  and  ^-direction  while  zero  forces  are  prescribed  in  the 
^/-direction;  these  boundary  conditions  enforce  symmetry  across  the  planes. 
The  atoms  in  the  y  =  1000  plane  (excluding  the  atoms  under  the  indenter) 
are  prescribed  zero  forces.  The  indenter  is  moved  downward  by  a  succession 
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of  increments  61  =  0.2  A  so  that  the  boundary  conditions  for  the  atoms  i  just 
below  the  indenter  are  given  by: 


Ui  =  (Q,-s6l,Q),  s  =  1, . . . ,  30  (18) 

Quasistatic  steps  are  then  considered  to  solve  for  the  displacements. 

The  site  energies  Ek{u)  of  each  atom  k  of  the  aluminum  crystal  are  modeled  by 
the  Embedded  Atom  Method  (EAM),  see  e.g.  [4,6].  Briefly,  the  semi-empirical 
potential  energy  for  atom  k  is  given  by 

Ek{u)  =  Fk(pk )  -  Ek\u),  (19) 

where  Fk  is  interpreted  as  an  electron-density  dependent  embedding  energy, 
pk  is  an  averaged  electron  density  at  the  position  of  atom  k,  and 

4V)  =  1  £  (20) 

(2) 

Here  Vk]  is  a  pairwise  potential  between  atoms  k  and  j  and  rkj  denotes  the 
interatomic  distance 

T'kj  ^/((*hfc  Xj)  -j-  (ltfc  H'j))  '  ((*®fc  T  ~{Uk  H'j'))  (21) 

Remark  4  ( Cutoff  functions )  The  quasicontinuum  method  employs  a  cut- 
off  function  to  approximate  the  interatomic  potentials  Vfd.  Since  the  potential 
decays  rapidly  with  respect  to  the  interatomic  distance  rkj,  the  potential  only 
includes  atoms  that  lie  within  some  short  distance  between  each  other. 

The  interatomic  distances  in  the  undeformed  configuration  of  the  crystal  are 
2.33  A  in  the  x- direct  ion.  One  (1,1,0)  layer  of  the  him  contains  about  1.3 
million  atoms  [19].  Note  that  the  lattice  is  two-dimensional,  but  that  dis¬ 
placements  are  allowed  in  three  dimensions  (constrained  by  periodicity  in  the 
^-direction)  and  the  energy  is  calculated  based  upon  three-dimensional  dis¬ 
placements. 

Rather  than  solving  for  the  solution  of  the  full  base  problem  (7),  an  “overkill” 
solution  of  the  surrogate  problem  (16)  is  considered  as  the  reference  solution. 
This  solution  hereafter  will  be  referred  to  as  the  base  model  solution.  This 
base  model  solution  involves  a  sufficiently  high  number  of  degrees  of  freedom 
so  that  it  is  considered,  for  the  purposes  of  this  study,  a  highly  accurate 
approximation  of  u. 

The  refinement  tolerances  were  set  to  0.000075  for  the  base  model  solution  and 
to  0.075  for  the  QC  solution  (this  value  is  recommended  by  the  authors  [10]). 
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The  meshes  corresponding  to  these  solutions  are  shown  in  Fig.  2.  The  vertical 
displacements,  at  an  early  stage  of  the  indentation  process,  and  just  before 
dislocation  nucleation,  are  shown  in  Figs.  3  and  4,  respectively.  In  the  early 
stages,  it  appears  that  the  displacements  computed  by  the  QC  solution  com¬ 
pare  well  with  the  base  model  solution.  However,  at  the  dislocation  nucleation, 
the  latter  appears  to  be  softer,  allowing  the  slip  plane  to  move  more  quickly 
into  the  material. 

For  better  comparison,  we  show  in  Fig.  5  the  magnitude  of  the  force  exerted 
by  the  crystal  onto  the  indenter.  This  force  represents  a  quantity  of  physical 
interest  as  it  clearly  indicates  the  nucleation  of  the  dislocation.  Note  that  the 
base  model  solution  has  been  run  only  to  load  step  27  due  to  computational 
cost.  However,  the  QC  solution  has  been  displayed  for  the  full  30  steps.  The 
QC  solution  seems  to  be  stiffer,  causing  the  dislocations  to  nucleate  one  load 
step  sooner  than  the  base  model  solution  as  the  critical  force  is  reached  more 
quickly.  Thus,  although  the  displacements  seem  to  compare  well  in  the  linear 
region,  small  errors  accumulate  resulting  in  an  inaccurate  portrayal  of  the 
mechanics  of  the  lattice.  Our  goal  in  the  next  sections  will  be  to  establish 
estimates  of  the  error  in  the  quasicontinuum  solution  with  respect  to  this 
quantity  of  interest  and  to  control  that  error  via  an  adaptive  algorithm. 

3  Error  estimation 

3.1  Errors  and  quantity  of  interest.  The  errors  in  the  QC  solution  u0,  with 
respect  to  the  solution  of  the  base  problem,  arise  from  three  sources:  (i)  use 
of  an  iterative  method  to  solve  the  nonlinear  problem,  (ii)  reduction  of  the 
number  of  degrees  of  freedom  from  N  to  R,  and  (iii)  approximation  of  the 
total  potential  energy  by  E0,  as  defined  in  (12)  and  (13). 

The  error  due  to  the  nonlinear  solver  is  controlled  at  each  iteration  and  is 
assumed  to  be  negligible  compared  to  the  other  sources  of  error.  This  error 
is  sometimes  referred  to  as  the  solution  error.  The  second  type  of  error  is 
analogous  to  discretization  error  in  Galerkin  approximations  such  as  in  the 
finite  element  method.  Here  it  can  be  regarded  as  a  model  reduction  error. 
Finally,  the  last  source  induces  a  so-called  modeling  error  due  to  the  modeling 
of  the  energy  using  the  coupled  local/nonlocal  QCM.  In  this  work,  we  will 
not  differentiate  the  three  types  of  errors  and  will  provide  for  estimates  of  the 
total  error. 

The  next  issue  when  dealing  with  a  posteriori  error  estimation  is  the  selection 
of  the  error  measure.  Early  works  on  the  subject  have  concentrated  on  global 
norms  such  as  energy  norms.  More  recently,  methods  have  been  developed 


10 


>■  -500 


1000 


Fig.  2.  Finite  element  triangulation  for  the  base  model  solution  (top)  and  the  QC 
solution  (bottom)  at  load  step  15.  The  base  model  solution  has  25484  atoms  while 
the  QC  solution  has  445  atoms. 
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3.  The  base  model  solution 


(left)  and  the  QC  solution  (right)  at  load  step  15. 
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Fig.  4.  The  base  model  solution  (left)  and  the  QC  solution  (right)  at  the  beginning 
of  load  step  27  and  26,  respectively.  The  base  model  solution  has  40554  atoms  while 
the  QC  solution  has  492  atoms. 


Fig.  5.  The  force-displacement  curve  comparing  the  evolution  of  the  base  model 
solution  and  the  QC  solution. 


to  construct  error  estimates  with  respect  to  quantities  of  physical  interest.  A 
typical  quantity  of  interest  in  the  nanoindentation  problem  (see  Section  2.3) 
is  the  reaction  force  on  the  indenter.  We  choose  here  this  force  as  the  quantity 
of  interest  for  determining  error  estimates.  Let  the  atoms  in  contact  with  the 
lower  surface  of  the  indenter  be  numbered  from  1  to  M.  Then  the  force  can 
be  written  as: 

M 

QH  =  ~Y.fi- n  (22) 

2—1 

where  n  is  the  outward  unit  normal  vector  to  <9Q  below  the  indenter.  In  terms 
of  the  potential  energies,  we  have: 

M  dE- 

Q(u)  =  n  (23) 

i= 1  OUl 
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We  note  that  Q(u)  is  a  nonlinear  functional  defined  on  the  solution  space  V. 
We  also  assume  in  the  following  that  the  meshes  are  constructed  in  such  a 
way  that  the  M  atoms  under  the  indenter  are  representative  atoms,  and  all 
the  forces  are  computed  by  the  nonlocal  approach. 

The  objective  is  then  to  estimate  the  error  quantity 

£  =  Q{u )  -  Q(iru0)  (24) 

where  nu0  E  V  is  obtained  from  u0  by  (10).  To  that  end,  we  follow  the 
approach  described  in  [11]  and  present  the  dual  problem  of  the  base  model  in 
the  next  section. 


3.2  The  dual  problem  and  error  representation.  The  general  approach  to 
obtain  estimates  of  the  error  £  makes  use  of  the  solution  of  the  dual  problem 
associated  with  the  base  model  (7): 


Find  p  E  V0  such  that 

B'(u;v,p)  =  Q'(u;v),  WveV0 


(25) 


where  the  derivatives  are  defined  as 


B\u]  v,  p)  =  lim  -  [B(u  +  6v ;  p)  -  B(u ;  p)} 

v— U 

Q'(u;  v)  =  lim  1  [Q(u  +  6v)  -  Q(u)] 


(26) 


In  the  molecular  statics  case,  we  have 


Na  Na 


B\u;v,p)  =  Y,Y, 


V; 


Q\u ;  v)  = 


3= 1  *=1 

Vi  ■ 


N 


^  d2Ek 

du.jdui  U 


Pi 


Na 


3= 1 


d2Et  ' 

dujdui y 


(27) 


However,  since  the  dual  solution  depends  on  the  exact  solution  u  and  that 
the  above  problem  may  be  intractable  due  to  the  large  number  of  atoms,  we 
may  use  instead  approximations  of  p.  An  approximation  can  be  obtained  by 
solving  the  surrogate  dual  problem: 


Find  p0  E  Wq  such  that 

B'0(u0;v,p0)  =  Q'0(u0-,v),  \/v  E  W0 


(28) 
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(29) 


and  by  extending  p0  to  the  space  V0  to  obtain  the  vector  7r 0p0  G  V0, 

R 

^oPo.i  =  J2Po,rM^i),  i  —  1,  .  .  .  ,  N 

r— 1 

The  quantity  of  interest  Q o  in  (28)  is  defined  for  v  G  W  by  Q o(u)  =  Q(ttv). 
In  the  following,  we  denote  by  e0  and  s o  the  errors  in  ttu0  and  7r0p0,  i.e. 


e0  =  u  —  7tu0  G  V 
S0=p~  7TqPq  G  V0 


(30) 


We  apply  the  results  of  [11]  to  derive  the  following  theorem  that  provides  for 
a  representation  of  the  error  in  the  quantity  of  interest: 

Theorem  1  Let  the  semilinear  form  £>(•;  •)  in  (7)  belong  to  CV{V)  and  let  the 
quantity  of  interest  Q(-)  as  defined  in  (23)  be  in  C3(0).  Let  u  e  V  andp  G  V0 
be  solutions  of  the  base  problems  (7)  and  (25),  respectively.  Let  (u0,p0)  G 
W  x  Wo  be  the  solution  pair  of  the  surrogate  problems  and  let  (7 m0,  rz 0p0) 
denote  their  extensions  to  the  spaces  V  x  Vo.  Then  the  error  in  Q(u )  produced 
by  7 Tu0  is  given  by 

S  =  Q(u )  -  Q(nuo)  =  TZ(nu0-,  p)  +  A(t tu0,  tt 0p0,  e0,  £0)  (31) 

where  TZ(nuo',v)  is  the  residual  functional, 

1Z{tzuo]  v)  =  F(v)  —  B(7vu0]  v),  veV0  (32) 

and  A  =  A  (71-1x0 .  e0i  £o)  Is  the  remainder, 

1  rl 

A  =  -  /  {B'fnuo  +  se0;  e0,  e0,  n  0p0  +  s£0) 

2  Jo 

—  Q"{tzuo  +  se0;  e0,  e0)}ds 

1  rl  (33) 

+  ~  /  {Q"'{ttuo  +  se0;  e0,  e0)  —  SB'finuo  +  <se0;  e0,  e0,  £q) 

2  Jo 

-  B"\nuo  +  se0;  e0,  e0,  e0,  7r0f»0  +  s£o)}(s  -  l)sds 


Goal-oriented  error  estimators  aim  at  estimating  £  by  accurately  approximat¬ 
ing  the  quantity  IZfirzuo',  p)  and  neglecting  the  higher-order  terms  A.  One  such 
approach  is  proposed  in  the  next  section. 

3.3  The  error  estimator.  We  first  rewrite  the  quantity  TMyiruo]  p)  in  different 
forms  in  order  to  lay  down  our  motivations  for  the  derivation  of  the  error 
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estimator.  Starting  from  the  definition  of  the  residual,  it  is  clear  that: 


TZ(nu0]  p )  =  F(ttu0,  p )  -  B(nu0;  p) 


Na 


Na 


=  E  fi  ■  Pi  -  E 


i— 1 
Na  ( 

=  E  (/< 

*=i  \ 

Na 


i=  1 
N 


N 


s-  QEk ,  v 

E  ■^r(7r«o) 


,fc=i 


du; 


Pi 


sr  dEk  <  \ 

E  ■^r(7rwo) 


fc=i 


duj 


Pi 


=  Eri(7mo)  'Pi 


2=1 


(34) 


where  the  residual  vector  r(nuo)  E  Vo  indicates  how  the  forces  acting  on  each 
atom  i  fail  to  be  equilibrated.  We  observe  that  the  calculation  of  Tl(7ru0-,p) 
may  be  cost  prohibitive  when  the  number  of  atoms  N,  or  rather  Na,  is  large.  In 
an  effort  to  reduce  the  computational  cost  of  the  error  estimator,  it  is  desirable 
to  take  into  account  only  those  contributions  that  are  the  most  significant, 
meaning  that  the  number  of  atoms  to  be  considered  for  the  calculation  of  (34) 
should  range  from  Ra  to  Na. 


Moreover,  thanks  to  the  linearity  of  the  residual  functional,  the  quantity 
lZ(iru0',p)  can  be  decomposed  as: 

TZ{nu0;  p)  =  TZ(ttu 0;  tt 0p0)  +  K(i ru0]  e0)  (35) 

It  is  well  known  that  the  contribution  1Z{tvuq\  7To p0)  vanishes  for  Galerkin 
approximations.  In  a  similar  manner  here,  this  term  fails  to  detect  the  model 
reduction  error.  It  follows  that  the  solution  p0  provides  for  a  poor  approxi¬ 
mation  of  the  dual  solution  p,  in  the  sense  that  TZ(ttuq;  tvoPo)  approximates 
72.(71 tuq,p)  poorly,  and  a  better  approximation  should  be  obtained  in  a  space 
larger  than  IV  q  • 


We  propose  here  to  evaluate  the  residual  and  the  dual  solution  on  a  mesh 
which  is  finer  than  the  mesh  used  for  the  evaluation  of  uq,  but  much  coarser 
than  the  mesh  that  would  be  obtained  by  considering  all  atoms  as  representive 
atoms.  Let  7 \  denote  such  a  partition  of  (we  shall  explain  below  how  to 
construct  7 \)  and  suppose  that  it  contains  a  total  number  of  N  nodes  with  Na 
interior  nodes.  We  introduce  the  vector  spaces  V  =  (Ma!)Ar  and  Vq  =  (WJ)Na 
as  well  as  the  extension  operators  n  :  V  — >  V  and  7T0  :  Vo  — >  Vo-  We  can  now 
define  a  residual  functional  1Z  on  V  such  that  for  any  u  E  V  and  v  E  Vq 

Na 

ll(u;v)  =  '52ri(u)  ■  Vi  (36) 

1=1 

where  the  r^u)  are  computed  via  the  coupled  local/nonlocal  quasicontinuum 
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approach.  We  will  also  consider  the  approximate  dual  problem: 

Find  p  e  Vo  such  that 

B'(nu o;  v,  p)  =  Q'(ttu0;  i),  Vi  £  Vo 


(37) 


with,  for  all  ii  G  V  and  v  G  Vo, 


Na 


B{u;v)  =  Y 


i=  1 


N  d  E, 


Y  Uk 

k=  1 


du 


[U] 


■  Vi 


~  "  dE i 

Q{u)  =  “  S  d^{u) ' n 


(38) 


Note  that  B  is  computed  using  the  N  representative  atoms  in  the  partition 
Vh-  We  emphasize  here  that  the  approximation  p  of  p  involves  the  same  types 
of  errors  as  in  u0,  but  that  it  also  strongly  depends  on  the  accuracy  of  the 
computed  solution  u0. 


We  now  define  the  error  estimator  with  respect  to  the  quantity  of  interest,  Q, 
as  the  computable  quantity 


Na 

V  =  Y^u0',  p)  =  Y  rii™ o)  •  Pi  (39) 

i— 1 

and  we  show  in  the  next  section  that  rj  is  a  reasonable  estimate  of  the  error 
S  =  Q(u)  —  Q 0(1*0)  =  Q(u)  —  Q{ttuq).  Finally,  the  finite  element  partition 
Vh,  or  enriched  mesh,  for  the  calculation  of  the  dual  approximation  p  and  the 
residual  r  is  constructed  using  the  adaptive  technique  described  in  Remark  3 
using  a  smaller  tolerance  than  0.075.  In  order  to  assess  the  quality  of  the  error 
estimate,  we  will  use  the  effectivity  index  defined  as  the  ratio  £  =  r]/£. 

3-4  Numerical  experiments.  We  perform  here  a  few  numerical  experiments 
using  the  same  setting  as  in  Section  2.3  in  order  to  study  the  performance  of  the 
error  estimator.  We  first  investigate  the  influence  of  approximating  the  dual 
solution  in  the  enriched  space  V  rather  than  in  the  space  Vq.  Close-up  views 
of  the  meshes  (QCM  and  enriched  meshes)  corresponding  to  these  spaces  are 
shown  in  Fig.  6.  As  expected,  the  error  estimator  performs  very  poorly  when 
the  dual  solution  is  approximated  on  the  QCM  mesh.  This  is  clearly  indicated 
in  Fig.  7  where  it  is  shown  that  the  error  estimator  detects  very  little  error  at 
all  load  steps.  By  contrast,  the  error  estimator  provides  reasonable  estimates 
when  the  enriched  space  V  is  used  for  the  approximation  of  the  dual  solution, 
as  shown  in  Fig.  8.  The  effectivity  indices  remain  mostly  close  to  unity,  except 
maybe  in  the  region  of  dislocation  nucleation  where  strong  nonlinear  behavior 
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Fig.  6.  QCM  mesh  (left)  and  enriched  mesh  (right)  at  load  step  9.  The  QCM  mesh 
has  432  atoms  while  the  enriched  mesh  contains  10887. 


Fig.  7.  The  (relative)  exact  error  and  the  error  estimate  on  the  QCM  mesh  without 
enrichment. 


occurs.  Indeed,  recall  that  the  QC  solution  dislocates  one  load  step  early. 
In  other  words,  the  primal  solution  uq  contains  large  errors  that  certainly 
pollute  the  approximation  of  the  dual  solution  at  that  particular  load  step. 
We  actually  show  in  Fig.  9  the  dual  solutions  p  and  p0  computed  using  the 
base  model  and  QCM,  respectively,  and  p  exhibits  many  more  details  than 
p0,  notably  in  the  region  away  from  the  indenter  near  the  slip  plane. 

4  Adaptivity 

We  propose  here  a  simple  adaptive  strategy  to  control  the  error  in  the  quantity 
of  interest  within  some  prescribed  tolerance  8toi-  Our  approach  is  different  from 
the  one  used  in  the  QCM  code,  but  was  made  to  fit  the  data  structure  available 
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Fig.  8.  The  relative  error  (left)  and  the  effectivity  indices  (right)  are  shown  com¬ 
paring  the  error  estimator  to  the  exact  error  for  the  enriched  QCM  mesh. 


Fig.  9.  The  dual  solution  of  the  base  model  (left)  and  the  QCM  (right)  at  the 
beginning  of  load  step  27. 


in  the  code. 

4-1  Adaptive  strategy.  For  the  purpose  of  spatial  adaptation  with  respect  to 
the  quantity  of  interest,  it  is  first  necessary  to  decompose  the  error  estimate 
77  into  local  contributions  that  could  be  employed  for  the  development  of  re¬ 
finement  indicators.  Due  to  the  structure  of  the  code,  we  have  adopted  an 
approach  in  which  the  local  contributions  are  defined  per  element  such  that: 

Ne 

V  =  dK  (40) 

K=  1 

This  is  accomplished  in  practice  as  follows:  for  each  element  K  of  the  partition, 
one  computes  the  nodal  contributions  rjf  =  f*j(7rito)-Pj  of  the  quantity  77,  from 
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which  one  can  calculate  an  elementwise  contribution  as: 


\I<\N, 


K 


.  Nd 

,  E  ‘f 

1k  i= 1 


c  f>i  d  x, 


(41) 


where  \K\  denotes  the  area  of  element  K ,  Nd  the  number  of  nodes  in  K ,  iV/4 
the  number  of  elements  sharing  node  i,  and  (f)l  the  restrictions  to  element  K  of 
the  linear  base  functions  as  defined  in  (9).  This  decomposition  is  simple  and 
easily  implemented  in  the  QCM  code,  but  is  not  unique. 


The  adaptive  algorithm,  the  so-called  Goals  algorithm  [13],  proceeds  as  follows: 

(1)  Initialize  the  load  step  to  s  —  0.  Input  user-tolerance  5toi- 

(2)  Go  to  the  next  load  step,  s  —  s  +  1. 

(3)  Solve  the  primal  and  dual  problems  as  described  in  Sections  2.2  and  3.2, 
respectively. 

(4)  Compute  the  error  estimate  as  discussed  in  Section  3.3. 

(5)  If  |)f|  <  fitoi\Q(uo)\,  go  to  step  (2).  Otherwise,  mark  those  elements  that 
satisfy  \t]K(ttu0,p)\  >  7 max  [i]K(nu0,p)\,  where  7  is  a  user-supplied 
number  between  0  and  1. 

(6)  Refine  flagged  elements  and  go  to  step  (3). 

Note  that  our  adaptive  algorithm  slightly  differs  from  QCM  in  the  sense  that, 
in  the  QCM,  the  elements  are  flagged  for  refinement  if  the  elemental  contribu¬ 
tions  are  below  some  user-specified  number  7 qc  and  that  the  adaptive  process 
within  each  load  step  eventually  ends  when  no  more  elements  are  flagged  for 
refinement. 


4-2  Numerical  examples.  In  the  following  examples,  we  choose  5toi  =  0.05 
(the  solution  is  controlled  so  that  the  relative  error  is  always  less  than  five 
percent)  and  7  =  0.25.  In  Fig.  10,  we  show  the  adapted  meshes  obtained 
using  the  QCM  and  the  Goals  algorithm  once  dislocations  have  nucleated.  As 
can  be  seen,  the  Goals  mesh  includes  many  more  atoms  near  the  indenter. 
Fig.  11  shows  the  evolution  of  the  number  of  atoms  (degrees  of  freedom)  for 
both  methods.  It  is  interesting  to  see  that  the  Goals  algorithm  adds  many 
repatoms  at  the  beginning  of  the  simulation  while  QCM  essentially  refines  at 
the  dislocation  nucleation. 

Force-displacement  curves  are  shown  in  Fig.  12.  We  observe  that  the  Goals 
algorithm  is  able  to  control  the  error  within  the  specified  tolerance  and  pro¬ 
vides  a  solution  that  predicts  the  dislocation  nucleation  just  as  the  base  model 
solution  does,  but  at  a  much  lower  computational  cost.  Relative  errors  and 
effectivity  indices  are  plotted  in  Fig.  13,  demonstrating  the  effectiveness  of 
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Fig.  10.  QCM  (left)  and  Goals  (right)  mesh  at  load  step  27.  The  number  of  atoms 
in  the  QCM  and  Goals  meshes  are  1629  and  3452,  respectively. 


Fig.  11.  Comparison  of  the  evolution  of  QCM  mesh  and  Goals  mesh. 


our  error  estimator.  Finally,  we  show  in  Fig.  14  the  dual  solutions  obtained 
using  the  base  model,  the  QCM,  and  the  Goals  algorithm.  Again,  this  result 
shows  that  the  dual  solutions  for  the  base  model  and  the  Goals  algorithm  are 
virtually  indistinguishable  while  the  solution  of  the  QCM  is  very  different. 

5  Conclusions 

In  the  present  work,  we  have  extended  the  methodology  of  goal-oriented  error 
estimation  and  adaptivity  to  molecular  statics  using  approximate  solutions 
produced  by  the  quasicontinuum  method  (QCM).  Estimates  of  the  error  in 
the  QC  approximations  with  respect  to  quantities  of  interest  are  derived  and 
are  used  as  a  basis  for  the  development  of  a  Goals  algorithm.  The  theoretical 
results  were  applied  to  a  sample  nanoindentation  problem  and  it  was  found 
that  the  Goals  methodology  provides  reliable  error  estimates  and  successfully 
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Fig.  12.  Force-displacement  curves  computed  from  the  base  model  solution,  the  QC 
solution,  and  the  Goals  solution. 


Fig.  13.  Relative  error  (left)  and  effectivity  indices  (right). 


controls  the  prediction  of  the  force  acting  on  the  indenter.  The  results  were 
consistently  verified  using  a  highly  resolved  solution  of  the  nanoindentation 
problem. 
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Fig.  14.  Comparison  of  the  dual  solution  for  the  base  model  (top  left),  the  Goals 
algorithm  (top  right),  and  the  QCM  (bottom)  at  the  end  of  load  step  27. 
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